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By means of molecular dynamics simulations based on the Billeter et al. [S. R. Billeter, A. 
Curioni, D. Fischer, and W. Andreoni, Phys. Rev. B 73, 155329] environment-dependent classical 
force field we studied the structural features of SiN^ samples at various stoichiometrics. Our results 
are in good agreement with experimental data and are able to reproduce some features which so far 
were not reproduced by simulations. In particular, we identified units containing N-N bonds, which 
are thought to be responsible for an unassigned peak in the radial distribution function obtained 
from neutron diffraction data and signals observed in electron spin resonance, X-ray photoemission 
spectroscopy, electron-energy-loss spectroscopy and optical absorption experiments. 

We have identified defects which are thought to be the responsible for the high concentration 
of charge traps that makes this material suitable for building non-volatile memory devices. We 
analyzed the dependency of the concentration of these defects with the stoichiometry of the sample. 
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I. INTRODUCTION 

Amorphous Silicon Nitride (a-SiaN4, or a-SiNx if the alloy is non-stoichiometric) is a very promising materials for 
application in the field of memory devices^ and optoelectronics^^. As for optoelectronics, in Si rich a-SiN^ samples Si 
clusters (Si-c) present in the matrbP are thought to be responsible for the luminescence of this materialP^l However, 
a-SiNz films are known to be luminescent matrices on their own due to the light emission related to the presence of 
specific defectsPtSl Si anc j N related defects act also as charge traps and their high concentration in a-SiN^ makes this 
material suitable for building non-volatile memory devices. 

Driven by its applicative importance, a significant effort has been made to clarify the atomistic and electronic 
structure of SiN^., and the relation between themP^On the one hand, the atomistic structure has been investigated 
by X-rajHand neutrorP diffraction experiments. These experiments revealed a significant degree of under-coordination 
of Si and N atoms, reporting and average coordination of 3.7 and 2.78 for Si and N, respectively. Another interesting 
structural feature found in the neutron Radial Distribution Function (RDF) is a peak at about 1.3 A. Such a short 
distance is only compatible with the presence of covalent bonds between N atoms in the sample. These bonds cannot 
be due to N 2 molecules trapped in the sample as in this molecule the bond length is shorter (1.1 A). An hypothesis is 
that this peak is due to units in which one N atom is bonded to another N atom and to Si atoms. As only the total 
RDF was measured, and the atomistic model used for its interpretation did not include N-N bonds, this hypothesis 
cannot be confirmed. As a result, the structure of SiN-r reported in Refs. [8| and |9] looks, locally, very similar to the 
crystal one. The presence of N-N bonds in hydrogenated and non-hydrogenated SiN x samples is further supported by 
Electron Spin Resonance (ESR)p21, X-ray photoemission spectroscopy, electron-energy-loss spectroscopy and optical 
absorption 1 -! experiments. 

Several computational studies have been performed to identify the atomistic structure and to compute the corre- 
sponding electronic structure of SiN x . Some of them used very rough atomistic models, such as the Bethe lattice^MI, 
in which the atomistic structure is decided a priori. We shall not compare our results with these works as their goal 
is to determine the dependence of the electronic structure from the atomistic configuration, without considering the 
statistical relevance of the (predetermined) configuration used in the calculation. In others works, the a-SiN x sa mples 
were generated by using procedures based on the randomization of either the atomic positions or the bonds^ES 
or on the combination of randomly oriented subunits 20 . Sometimes the random-orientation step was followed by a 
relaxation s tep perform ed by low-temperature Molecular Dynamics (MD) or Monte Carlo (MC) runs. Finally, in 
some studies^ * 16 * 18 * 21 ! the samples were prepared by the quench- from-the-melt procedure (see Sec. [n]). While all 
these studies are able to reproduce the general features of the experimental RDF, in none of them it is reported 
the formation of N-N bonds. Only the pair distribution function reported in Ref. [T7] shows indeed a small peak 
at ~ 1.3 A, but this feature is not discussed by the authors. Moreover, its intensity is very low compared with 
experimental dataP The failure of classical and, especially, ab initio force models to reproduce this covalent bond is 
surprising. We believe that the reason of this failure has to be sought in the procedure by which the samples were 
generated rather than in the accuracy of the force models used for the simulations, with the exception of Ref.s [2] 
and [18] in which it is used a force model that explicitly prevents the formation of N-N. In fact, in the samples 
produced by random orientation of subunits the formation of the N-N bonds is prevented by the procedure that 
imposes t he con servation of the Si/N alternation. As for the samples prepared by ab initio quenching- from-the-melt 
procedure s] 15 * 16 ^, the melting temperature is too low (3000 — 3500 K) and the melting and cooling times are too short 
(few picoseconds each). In these conditions the breaking of Si-N bonds and the formation of N-N bonds, required 
for the formation of units involving N N bonds, cannot take place. In fact, for a similar material (a-SiOa) Vollmayr 
et al|22l have demonstrated that higher temperatures (7000 K) and very long melting and cooling processes (up to 
100 ps and 1 ns, respectively) are needed to obtain good amorphous samples. A special consideration is in order for 
Refs. [18] . In fact, the results reported in this paper are obtained using a potential in which the N-N interactions 
are purely repulsive and therefore the formation of N-N bonds is prevented. Also Vashishta and KaliaP^ did not find 
N-N bonds in their samples generated via classical MD simulations. In this case the melting temperature, the melting 
time and the cooling rate are adequate. In this case, we believe that the inability to produce units containing the 
N-N bond is probably due to a limitation of the force field adopted for the simulations. 

The methods and the conditions used for generating the samples could have introduced also other artifacts in the 
corresponding atomistic structures. For example, they could have produced samples with a too low concentration 
of Si and N dangling bonds, over-coordinated Si and Si-Si bonds, which are all thought to be defects influencing 
the electronic and optical properties of the material. The aim of this work is therefore to further investigate the 
atomistic structure of a-SiN K at various stoichiometries (0.3 < x < 1.8) using the quenching-from-the-melt procedure. 
We shall extensively explore the domain of the parameters governing this procedure: i) melting temperature, ii) 
duration of the melting and iii) quenching rate. In addition to a general agreement with the experimental RDF, 
we shall use the formation of N-N bonds, and correspondingly the agreement with neutron RDF at ~ 1.3 A, as a 
benchmark for assessing the quality of the samples obtained. We shall further validate these structures by performing 
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ab initio geometry optimization and MD simulations. On good samples we shall also identify the other relevant 
structural parameters, such as concentration of coordination defects, and concentration and structural parameters of 
the N(Si m N n ) and Si(Si m N„) subunits contained in them. 

In our stu dy we perform classical MD simulations based on the environment dependent force field recently developed 
Billcter et al! 24 * 25 l The ability of this force field to model the equilibrium structure and the dynamical properties of 
the closely related a-SiC>2 system has already been established 2 ^ 2 ^. This makes us confident that this potential is 
also able to properly represent the interatomic interactions in SiN^ as well. However, to the best of our knowledge, 
no works have been published before which make use of this potential to study a-SiN x . Therefore, another objective 
of this research is the validation of this potential also for this material. 

The paper is organized as follows: in Sec. [n]we give the computational details of our simulations. In Sec. |III| we 
shall present the results of our simulations and the analysis of the atomistic structure of the a-SiN^ samples. Finally, 
in Sec. IIVI we shall draw some conclusions. 



II. COMPUTATIONAL DETAILS 



In this paper we generate and analyze the structure of amorphous SiN^ sample at v arious compositions as obtained 
from MD simulations based on the Billeter et al. environment dependent force field 124125] This force field allows to 
treat silicon oxide and oxynitride systems, both hydrogenated and non hydrogenated. The Billeter et al. force field 
is an evolution of the Tersoff potential and one of its advantages is that it describes more accurately coordination 
defects. In fact, one problem with the original Tersoff potential was the tendency to produce structures with an 
high abundance of coordination defects. In the Billeter et al. potential this problem is solved by adding a "penalty" 
term that depends on the deviation from the optimal coordination. Moreover, with specific reference to the case of 
aSiN x and at a variance form previous classical force field (e.g. the de Brito Mota et al. force fielcP^]), the Billeter 
et al. potential takes explicitly into account the formation of N-N bonds. The reliability of this force field to model 
complex systems and phenomena, such as Si/SiC>2 interface^, Si nanoparticles embedded in amorphous SiOjf^, Si 
and O self-diffusion in SiO^^ and strain distribution in Si nanowires^ has been already proven. On the contrary, to 
the best of our knowledge, a careful validation on realistic silicon nitride systems is still lacking. This validation will 
therefore be another objective of this work. 

The a-SiaN4 sample was obtained by means of the quenching-from-the-melt method. This technique consists 
in melting a crystal of proper stoichiometry at high temperature and then slowly cooling it down to the target 
temperature, at which all the relevant properties are computed. The cooling step of the process can be performed in 
different ways, for example by a constant temperature MD in which the target temperature is continuously decreased or 
by running a series of consecutive constant temperature MD simulations in which the target temperature is decreased 
when passing from a run to the next. In this paper we have used the first approach (continuous cooling). However, we 
have also prepared a selected set of samples following the stepwise cooling procedure. The comparison of the structural 
characteristic of the samples prepared with the two methods has shown that there is no significant difference between 
them provided that the average cooling rate (i.e. dT/dt) is the same. Summarizing, the parameters controlling the 
quenching-from-the-melt procedure are: the melting temperature T m , the duration of the melting step t m and the 
quenching rate dT/dt. The effect of these parameters on the structural properties of the amorphous system has been 
evaluated by producing several samples according to the values reported in Table [TJ These results are described in 
detail in the next section. 

The initial configuration for the quenching-from-the-melt procedure for the stoichiometric sample was the crystal 
structure of the a-Si 3 N 4 phase. The density was kept fix to the experimental density of the a-Si 3 N 4 (p = 3.14 g/cm 3 ^, 
which is slightly lower than the density of a-Si3N4 (p ~ 3.2g/cm 3 ). In order to check for finite size effects, we ran 
simulations on samples containing as many as 224, 1792 and 14336 atoms without observing any significant difference 
between the structure (RDF, angular distribution function and coordination defects) of the samples. This means 
that there is no correlation beyond the maximum length-scale compatible with the simulation box of the 224 atoms 
sample, that is ~ 5. 8 A. This is indeed in agreement with the observation that the g(r) of th e am orphous samples do 



not show any correlation beyond 3.5 A (see next section). All the results reported in Sec. Ill refer to the smallest 
sample. 

The amorphous sub(supra)-stoichiometric samples were obtained starting from the a-SisN4 sample and randomly 
replacing N (Si) atoms with Si (N) atoms. In particular, we generated SiN x samples at x = 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.8. 
The simulation box was isotropically modified so as to obtain a density consistent with the experimental density,^ 
which varies with the composition (see Table [II]). The so obtained samples are melt at T m = 6000 K for t m — 100 ps 
and then quenched down to room temperature with a quenching rate of dT/dt = 5 K/ps, that is at the lowest 
quenching rate used for generating the stoichiometric sample. 

We further checked the reliability of our results by comparing the theoretical density of all the samples at ambient 
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conditions with experimental data.^ To this end, we performed constant pressure MD simulations according to the 
Martyna- Tobias-Klein methocP^ for non-isotropic cells starting from the equilibrated samples at fixed volume (see 
above). 

Finally, on samples at composition x = 0.4, 1.333 and 1.8, i.e. in the sub-stoichiometrical, stoichiometrical and supra- 
stoichiometrical domains, we performed ab-initio geometry optimization and MD simulations. These calculations have 
been performed using a setup as close as possible to the one used by Billeter et al. to generate the training set used 
for fitting the parameters of their force field. In particular, the calculations have been performed using the version 
3.13.2 of the CPMD simulation packaged The electronic structure calculations were performed using the Density 
Functional Theory in the generalized gradient approximation of Perdew-Burke-Ernzerhof. 34 The Kohn-Sham orbitals 
were expanded on a plane wave basis set with a 80 Ry cutoff. Finally, electron-ion interactions are described using 
Troullier- Martins pseudopotentialsP^ 



III. RESULTS AND DISCUSSION 

The general aspects of the quenching-from-the-melt amorphization protocol were presented in sec. [TTJ. The ability 
of the protocol to produce a realistic amorphous sample depends on the parameters controlling the various steps. In 
particular, the parameters must be chosen such that in the first step the sample is completely melt and in the second 
step the sample can relax to a local minimum of the free energy. In this work we have ran a series of simulations 
at various T m , t m and dT/dt. Moreover, we have evaluated the effect of the stoichiometry on the structure of the 
sample by studying, in addition to the stoichiometric SiaN4 sample, sub and supra stoichiometric SiNj, samples, with 
x = 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.8. In table[l]we summarize the list of simulations performed on the SiaN4 sample and the 
corresponding values of the parameters governing the quench-from-the-melt procedure. 

Before presenting and discussing our results, it is worth to comment on the melting temperatures used in our 
simulations. The experimental melting temperature of the Si3N 4 depends on the conditions under which the sample 
is prepared and the melting temperature is measured. Commercial samples (e.g. a-Si3N 4 from Atlantic Equipment 
Engineers, Div. of Micron Metals, Inc. - Bergenfield, NJ) are reported to melt at about 2175 K (other sources report 
a-SisN4 to decompose before melting). Accurate experiments^ in N2 atmosphere have shown the a-Si3N4 to melt 
approximately at 2775 K. However, in simulations the melting process typically occurs at much higher temperatures. 
There are several reasons for this, of which the accuracy of the force field is possibly not the most relevant one. In 
simulations we typically start from the perfect crystal while in real samples there is a natural abundance of defects, 
such as vacancies, which allow a larger atomic mobility. Moreover, especially in constant volume MD simulations 
of small samples, the fluctuation of local density is strongly limited. Fluctuations of the density naturally occurring 
in real samples enhances atomic diffusion which, in turn, produces the melting. The higher temperature needed for 
observing the melting in simulations might also be due to another phenomenon, namely that the formation of a liquid 
in a bulk crystal (i.e. in absence of a surface, as it is the case in all the simulations mentioned in the introduction and 
in the present work) this process is termally activated (i.e. there is a barrier in the free energy profile separating the 
crystal from the liquid domain) . An exhaustive description of the theory of the nucleation of a liquid phase in a bulk 
crystal, and the relative experimental and computational results, is given in the Chapter 3 of Ref. 1371 This means 
that at the computational melting temperature (i.e. the temperature at which the liquid phase is more stable than 
the crystal phase given the specific force field) the system can stay trapped in a metastable state corresponding to the 
overheated crystal for the entire duration of the simulation. Therefore, in order to obtain a liquid sample starting from 
a crystal on the short timescales allowed by simulations (from few ps in ab-initio simulation up to several hundred ps 
in classical MD with environment dependent classical force fields) it is necessary to go much beyond the experimental 
or computational melting temperature. This explains why in our simulations we perform the melting at very high 
temperatures (all above 5000 K). This is in agreement with what was done in analogous studies on a related material 
using a different force fields 

In Fig. 1 it is reported the pair correlation function (<?(r)) of the SisN4 melt at various temperatures measured 
over 10 ps of constant temperature MD after initial 100 ps of thermalization at the target temperature T m . For 
comparison, we also report g(r) of the SisN4 crystal at 300 K. In the inset of Fig. 1 we report the N-N partial g(r) 
at the same temperatures. The first remark is that for T m < 5500 K the structure preserves a long range order. 
For example, these samples still have peaks in the g(r) at ~ 3.4 A and ~ 4.1 A and a series of superimposed peaks 
between ~ 4.3 — 4.8 A. Moreover, the broad peak at ~ 2.8 A still shows traces of the three peaks present in the crystal 
structure. This means that temperatures < 5500 K do not produce adequate melting on a timescale of 100 ps. In 
these samples, there is no evidence of the peak at ~ 1.3 A. However, for slightly higher temperatures (T m > 5800 K), 
the long range features in the g(r) disappear. This means that the sample is completely melted on the 100 ps timescale 
of the simulation. In these liquid samples it also appears the peak at 1.3 A. The shape of the g(r) in general, and the 
peak at 1.3 A in particular, are only marginally affected by a further increase of the temperature. This indicates that 
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the liquid samples produced by high temperature simulations have reached the equilibrium. 

We now turn to the analysis of the origin of the peak at 1.3 A. Such a short distance is compatible only with N-N 
bonds. This is confirmed by the presence in the N-N partial g(r) of a peak at 1.3 A (see inset in Fig. 1). This result 
is consistent with experimental dataP^ The presence of N-N bonds indicates that in the melt the N atoms are not 
only part of the N(Si 3 )/Si(N 4 ) units, which are the only units present in the crystal. Let us indicate by N(Si„N m ) a 
generic unit in which m nitrogen and n silicon atoms are bonded to a "central" N atom, and by Si(Si n N m ) a unit in 
which m nitrogen and n silicon atoms are bonded to a "central" Si atom. The concentration of the N(Si„N m ) units 
can be used to monitor the stoichiometric defects of the central N atom and the concentration of the Si(Si n N m ) units 
those of the central Si atoms. In table |Ln| is reported the concentration and the structure of the more abundant of 
such units at 5800 K, 6000 K and 7000 K. Let us start by analyzing the stoichiometric defects of the N atoms. In 
the melt at 5800 K, the most abundant N(Si3)/Si(N4) units are (in descending order of concentration) N(Sia) (58%), 
N(Si 2 N) (23%), N(Si 2 ) (13%), N(SiN) (3.5%) and N(SiN 2 ) (3.5%). Consistently with what we found for the g(r), 
the temperature affect only marginally the abundance of these units. At 6000 K we found a concentration of 48%, 
23%, 14.5%, 4.5% and 4% for the units, N(Si 3 ), N(Si 2 N), N(Si 2 ), N(SiN) and N(SiN 2 ) respectively. At 7000 K the 
concentration of these units are (in the same order as before) 41%, 22%, 18.5%, 6.5% and 4.5%. As expected, the main 
effect of the temperature is to increase the concentration of defective units (N(Si 2 N), N(Si 2 ), N(SiN) and N(SiN 2 )), 
some of which contain N-N bonds, and to decrease the concentration of N(Sia), i.e. the unit present in the crystalline 
SiaN4. As for the structure of these units, the Si-N bond length ranges from 1.77 A to 1.90 A, and the N-N one 
from 1.41 A to 1.53 A depending on the unit. As expected, the bond lengths is shorter in under-coordinated units 
(n + m < 3) and longer in over-coordinated ones (n + m > 3). On the contrary, the stoichiometry of the structure 
affects only marginally the structure. The situation is very similar for the Si(Si„N m ) units: the concentration of 
defective units, some of which containing Si-Si bonds, increases with the temperature. Also in this case, the structure 



is only marginally affected by the temperature. These results were confirmed by ab initio simulations. In table III 
we compare the structural data and concentration of the N(Si„N m ) units as obtained from our classical and ab initio 
MD simulations. In particular, the latter were started from one configuration extracted from the classical simulations 
at corresponding temperatures. The comparison shows that both the concentration and the geometry of the classical 
and ab initio N(Si n N m units are very close to each other, so confirming the reliability of the a-SiN3 structure obtained 
from the classical simulations. 

Also the quenching rate is expected to influence the structure of the a-Si3N4, we therefore studied its effect by 
performing a series of simulations at various dT/dt in the range 0.5 K/ps to 500 K/ps (see Table [!]). For these 
simulations we started always from the melt generated at T m — 6000 K. In Fig. 2 are shown the g(r) obtained 
from simulations at various quenching rates. The main features of these g(r) are very similar, namely all the g(r) are 
characterized by three main peaks at 1.3 A, 1.7 A and 2.9 A, and another peak at 2.65 A merged with the previous one. 
However, there are significant differences in the details between the g(r) at the various quenching rates. The intensity of 
the N-N peak decreases and the intensity of the two main peaks at 1.7 A and 2.9 A increase for slower quenching rates. 
This dependence of the g(r) on the quenching rate can be quantified by computing the number of pairs corresponding 
to the different peaks. In turn, these can be computed by properly integrating the g(r) on the intervals associated 
to these peaks (i.e. taking into account the Jacobian associated to the cartesian-to-polar coordinate transformation). 
However, the peak at 2.9 A is the combination of three superimposed peaks associated to second Si-Si, Si-N and 
N-N neighbor pairs. This peak is therefore not suitable for the analysis proposed above. For this reason we perform 
the analysis described above by comparing only the number of N-N and Si-N pairs corresponding to the peaks at 
1.3A and 1.7A respectively. In the inset of Fig. 2 we show the ratio of the number of Si-N pairs over the number of 
N-N pairs as a function of the decreasing cooling rate. This curve clearly demonstrates that this ratio increases with 
the decreasing of the quenching rate, eventually converging to an equilibrium value for very slow quenching rates (for 
dT/dt < lOK/ps). This behavior can be explained as follows. Let us consider the N(Si 3 ) and the Si(N4) units, that 
are the building blocks composing the (stoichiometric) Si3N4 crystal. Assuming that these units are more stable than 
the units containing the N-N bond (N(Si„N m ), see Fig. 3), which is a reasonable assumption as the second is not 
present in crystalline (low temperature) Si3N4, the relative concentration of N(Si n N m ) units with respect to N(Si 3 ) 
will increase with the temperature. Therefore, at 6000 K there will be a higher relative concentration of these kind of 
units with respect to the equilibrium concentration in a-Si3N4 at 300 K. When the samples are quenched too quickly 
the system cannot relax to the equilibrium and the final relative concentration of N-N pairs is too high, i.e. it will 
remain close to the concentration present in the high temperature melt. On the contrary, at low quenching rates the 
system can relax and the N-N pairs concentration can reach its equilibrium value. In Fig. 4 we report the Radial 
Distribution Function (RDF) of the sample obtained at dT/dt = 5 K/ps, which is in very good agreement with the 
neutron diffraction data reported in Ref. [9j In particular, the agreement of both the position and the intensity of the 
computational and experimental peak at ~ 1.3 A is very remarkable. However, as a further validation of our results, 
starting from the last configuration of the classical MD simulation we performed an ab-initio geometry optimization 
of the structure. As a first remark, all the N-N bonds present in the initial configuration (i.e. the one obtained from 
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the classical MD) are still present in the final one. Moreover, the root mean square displacement (RMSD) between 
the ab-initio minimum energy configuration and the classical configuration averaged over the MD trajectory (there 
is no diffusion at this temperature over the timescale of the simulation and therefore this averaged configuration is 
meaningful) is ~ 0.1 A. This result indicates that the state identified by classical MD corresponds to a metastablc 
state also on the (more accurate) ab-initio potential, so confirming the reliability of our simulations. Finally, in order 
to take into account the anharmonicity of the potential, and the possibility that the metastable state identified by the 
classical simulation is not the most stable one, we ran a 5 ps long ab initio MD simulation. In Tables |IV| and [V] are 
reported the concentration and the structure of the N(Si„N m ) units as obtained from the classical and ab initio MD 
simulations, respectively. By comparing the results reported in these two tables, we notice that the main difference 
among the classical and ab initio samples is the concentration of the N(Si4) unit, which is lower than 1% in the former 
and about 7% in the latter. As for the bond length, the maximum difference between the two samples is ~ 6%. 

Analyzing more in detail about the structure of the N(Si„N m ) units, we found that the N-N and Si-N bond length 
are rather unaffected by the composition of the unit as far as the coordination number of the central N atom is 
preserved (i.e. m and n at a fixed n + m), being N-N and N-Si bond length ~ 1.35 A and ~ 1.75A, respectively. 
Also the bond angle is essentially unaffected by the chemical composition of the N(Si„N m ) units and its value is 
~ 120°. Such a bond angle indicates that the structure of these units is trigonal. These results are in agreement with 
experimental dataP In conclusion, in the stoichiometric sample the N-N bonds are single bonds occurring in trigonal 
units at various concentration of Si and N. This result give evidence that the assumption that if two N atoms bond 
together they will form a N2 molecule (that would leave the sample), which is the assumption made in the de Brito 
Mota et al. force fielcP^ to make vanishing the N-N attractive interaction, is wrong. Below we shall show that this 
conclusion is independent on the stoichiometry of the sample. 

Before moving to the analysis of samples at different stoichiometries, we still have to evaluate the effect of the 
quenching rate on the concentration of Si and N coordination defects. This is a further validation of our computational 
setup and it is of great importance as the concentration of these defects is considered the structural features at the 
base of the interesting optical and electrical properties of this material. We do this by computing the percentage of 
Si" and N m in the sample, with n = 3,4,5 and m = 2,3,4, where the superscript indicates the number of bonds 
formed by a given atom. The coordination number of a given atom is computed according to the distance between 
this atom and its neighbors: two atoms are considered bonded if their distance is lower than the first minimum of 
the corresponding partial g(r). In Fig. 5 it is shown the percentage of Si 11 and N m as a function of the cooling rate. 
These data show that the concentration of N m is only marginally affected by the quenching rate. Essentially, N atoms 
are almost all 3- fold coordinated at any quenching rate apart for the case of the very high quenching rate dT/dt = 
500 K/ps and dT/dt = 100 K/ps. The effect of the quenching rate on the coordination defects of Si is much more 
evident. At very high quenching rates the percentage of 5, 4 and 3-fold coordinated Si atoms is ~ 22 %, ~ 71 % and 
~ 7 %, respectively. The concentration of defects initially increases, then suddenly decreases with lower quenching 
rates and finally the sample converges to a composition of %Si 3 ~ 24 %, %Si 4 ~ 75 % and %Si ~ 1 %. It might 
appear surprising that the percentage at high and low quenching rates are close to each other. We believe that this 
is just due to chance. 

Based on the analysis of the effect of the melting temperature and quenching rate on the stoichiometric sample, 
the non-stoichiometric samples were prepared with a melting temperature T = 6000 K and a quenching rate dT/dt = 
5 K/ps, that guarantees converged g(r) and composition of coordination defects. As in the case of stoichiometric 
sample, we checked the reliability of our results by performing ab-initio geometry optimization and 5 ps long MD 
simulations. Due to the high computational cost of ab-initio calculation on samples containing more than 200 atoms, 
we limited this analysis to x = 0.4 and x = 1.8 samples only, i.e. one in the sub and one in the supra-stoichiometric 
domain. Similarly to what found for the stoichiometric sample, the (averaged) classical and ab-initio configurations 
resulted to be very close to each other for both the sub and supra-stoichiometric samples (RMSD ~ 0.2 in both cases). 
This means that the state found by classical MD simulation is a metastable state also on the ab initio potential 
energy surface. We also compared classical and ab initio results at finite temperature (see Tables IV and |v|. Both 
the concentration of the N(Si„N m ) units and their structure as obtained from classical and ab initio are in agreement. 
This confirm that the results obtained from the classical simulations are reliable also in non-stoichiometric conditions. 
Finally, in order to further compare our results with experimental data, we also performed NPT MD simulations at 
room temperature and ambient pressure aimed at computing the density of all the samples. These results are reported 
in Fig. 6 together with experimental data. The latter have been obtained by extrapolating experimental densities 
published in Ref. [llj to the stoichiometries studied in the present work. As a general remark, the computational 
trend is very similar to the experimental one, with the density growing with the concentration of N in the sample. 
More in detail, the computational values are very close to the experimental one, being the difference between ~ 1.5 % 
and r-j 5 % (see the inset of Fig. 6) . 

Moving to the analysis of the structure of non stoichiometric samples, Fig. 7 shows the total, and Si-Si and Si-N 
partial g(r) for SiN^, with x=0.4, 0.6, 0.8, 1.0, 1.2, 1.33 (S13N4), 1.4 and 1.8. The N-N partial g(r) is not shown as, 
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apart for the peak at 1.3A, its shape is very broad and no insight on the structure of the samples can be obtained 
from it. In the following we shall analyze the peaks at 1.3 A, 1.7 A and the complex group of peaks between 2.2 and 
3 A. As for the first peak, already identified as due to N-N pairs in the stoichiometric sample, it increases with the 
nitrogen content. While this trend is in general not surprising, it is interesting to notice that this peak is already 
present in sub-stoichiometric conditions, starting from samples of composition x = 1. The main peak at 1.7 A, due 
to Si-N pairs, initially increases with x reaching a maximum at x = 1.33 (stoichiometric composition) and then 
decreases. This is because at low x there are not enough N to match with Si and the situation is reversed for x > 1.33. 
More interesting is the evolution of the complex series of peaks between 2.2 and 3 A. The g(r) in this interval is 
due to the superposition of the contributions of the three pairs N-N, Si-Si and Si-N. However, at low x this peak is 
mainly due to Si-Si pairs, as it is confirmed by the comparison of the total and Si-Si partial g(r). With the increase 
of N content the shape of the g(r) in this range becomes more complex. At very high x a satellite peak at 2.7 A 
appears. This satellite peak is due to Si-N pairs, as can be seen by comparing the total and Si-N partial g(r). As for 
the dependency of the Si-Si partial g(r) with the composition, we notice that the peaks at 2.3 A, due to chemically 
bonded Si atoms, decreases in intensity with the increase of x. At x = 1.33 this peak is very low, even though it never 
vanishes completely even at very high x. On the contrary, the intensity of the peak at 2.9 A increases with a;. This 
latter peak is due to non bonded Si-Si pairs of atoms which are both connected to the same N (see top-left panel 
of Fig. 7). The opposite trend of the two peaks can be explained by the transformation of Si(Si m ) units, which are 
present in sub-stoichiometric samples, into N(Si m ) units, typically present in stoichiometric and supra-stoichiometric 
systems. Very interesting is the analysis of the origin and the trend of the peak at 2.6 A in the Si-N partial g(r). 
This peak is due to non-bonded Si-N pairs belonging to a unit in which the two atoms are both bonded to a N (see 
the pair connected by the dashed line in Fig. 3/B). It is worth to mention that this peak appears for the first time at 
x = 1.0, in correspondence with the appearance of the peak at 1.3 A, and then increases with x. On the basis of this 
analysis we conclude that in SiN^ samples with x > 1 are presents units of the kind N(Si„N m ), that is units in which 
a N is bonded to both Si and N (see Fig. 3/A-C). Similarly to the case of the stoichiometric sample, we analyzed 
the nature, concentration and structure of these units (see Table IV). In sub-stoichiometric samples we only found 
N(Si2N) units. This is, of course, due to the low concentration of N. On the contrary, in the stoichiometric and in 
supra-stoichiometric samples we found N(Si2N), N(SiN2) and N(Na) units (the latter was found only in the samples 
corresponding to x — 1.6 and x — 1.8). The concentration of N(Si ra N m ) units in general, and the concentration of 
to > 2 units in particular, increases with x. The total concentration of these defective units is already very significant 
in the stoichiometric sample (~ 16.5%) and becomes higher than the concentration of N(Si3) units at x = 1.8 (58% 
vs 39.5%). This latter result is confirmed by the ab initio MD simulations in which the concentration of N(Si„N m )) 
units resulted to be ~ 53%, to be compared with a concentration of the N(Sis) unit of ~ 37%. Moving to the analysis 
of the structure of the N(Si„N m ) units, both classical and ab initio simulations indicates that it is only marginally 
affected by the stoichiometry of the sample and the chemical composition of the unit. On the contrary, as expected, 
it strongly depends on the coordination of the central N atom. The average bond length in N(Si n N m ) units in which 
the central N is three-fold coordinated (i.e. n + m = 3) is cLnn — 1-33 A and d,NSi = 1-74 A, and d^N — 1.42 A 
and oInSi — 1-8 A in classical and ab initio simulations, respectively. As for the (average) bond angle (a), classical 
simulations give always an (almost perfect) trigonal structure for all the units in which the central N is three-fold 
coordinated (a ~ 119°). Ab initio simulation show a larger departure form this structure, with an (average) bond 
angle of ~ 116°. 



Another very relevant structural parameter to consider is the concentration of coordination defects. In fact, as 
already mentioned in the introduction, both Si and N coordination detects are thought to be responsible for the 
electronic properties of a-SiNj^Sl l n Fig. 8 we report the concentration of Si" (top) and N m (bottom) as a function 
of x. As for N, we did not observe any N 2 at any stoichiometry. On the contrary, there is a high concentration of N 4 
defect (thought to be an electron trapping defect) at low x. However, the concentration of this type of defect decreases 
with x and becomes negligible for x — 1.2. The situation is different for Si. At x = 0.4 the concentration of Si 3 , Si 4 
and Si 5 is 15 %, 79 % and 6 %, respectively. By comparison with the concentration of similar defects in amorphous 
Si (a-Si), we conclude that the effect of low concentration of N is to increase the number of under-coordinated Si. In 
fact, in a-Si, by using the same potential and amorphization protocol, we have measured 3% of Si 3 , 93% of Si 4 and 
4% of Si 5 (these values are in line with ab-initio MD simulation results 39 ). At variance with N, the concentration 
of the Si coordination defects remain large at any N content. However, starting from the stoichiometric composition 
we observe a reduction of their concentration. Much the same as for N, most of the Si coordination defects are due 
to under-coordination. However, at very low and very high N content there is also a non negligible amount of 5-fold 
coordinated Si atoms. 



IV. CONCLUSIONS 

We studied the structural features and the coordination defects of silicon nitride samples at various compositions 
as obtained fr om the qucnch-from-the-melt method. Our results are in good agreement with X-ray and neutron 
diffraction datcP^I. In particular, our simulations show that the peak at 1.3 A in the experimental g(r) is due to 
the formation of N(Si n N m ) units. This finding is consistent with other experimental data, namely Electron Spin 
Resonance (ESR) measures^, and X-ray photoemission spectroscopy, electron-energy-loss spectroscopy and optical 
absorption 11 experiments on hydrogenated and non-hydrogenated SiN x samples. Also the density at room temperature 
and ambient pressure obtained from NPT-MD simul ations is in good agreement with experimental results.^ This 
is a further confirm that the Billeter et al. force fiel d 24 * 25 * is adequate for simulating semiconductors of SiN^ type. 
We believe that the inability of previous computational works to reproduce the features mentioned above was due 
to either the force field used for the simulations, which explicitly prevented the formation of N-N bonds, or to the 
quenching-from-the-melt parameters used in the simulations. In particular, we think that in previous simulations^^ 
it was used a too low melting temperatures and/or a too short melting time. Since the breaking (of Si-N) and the 
formation of (N-N) bonds is a rare events, i.e. it occurs on a time scale longer than the typical duration of an MD 
simulation, the conditions used in previous simulations do not allow them to occur on the time scale of simulations. 
Other structural features, and as a consequence electronic structural features, might be also affected by the inadequacy 
of the force field or the amorphization protocol. 

We also analyzed the effect of the quenching rate on the structure of the sample, concluding that very low quenching 
rates are needed to obtain well converged amorphous samples. 

We also studied the structural features of non-stoichiometric SiN x samples, which are of great technological interest. 
We found that at low x the sample contains Si(N m Si4_ m ) and NSi3 units. However, already at x = 1, i.e. still in 
sub-stoichiometric conditions, other units containing N-N bonds are formed. These units are not N2 molecules, rather 
N(Si n N m ) units in which a nitrogen atom is bonded to both Si and N atoms. 

Finally, we studied the concentration of coordination defects as a function of the stoichiometry of the sample. These 
defects are considered relevant for the electronic properties of this material, especially for its use in the non-volatile 
memory device field. We found that at very low N concentration there is a significant amount of 4-fold coordinated 
N atoms while the concentration of N 2 defects is negligible at all x. Defects of this type are thought to be responsible 
for the electron trapping properties of SiN x . As for Si, we found that there is a significant amount of coordination 
defects at any N content. In particular, these defects are mostly associated to under-coordinated Si, apart at very 
low and very high x where there is also a non negligible amount of Si 5 . 
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FIGURE CAPTIONS 

Fig. 1: Pair correlation function of liquid stoichiometric samples at various melting temperatures after 100 ps of 
simulation. In the inset are displayed the partial N-N g(r) at the corresponding melting temperatures. 

Fig. 2: Pair correlation function of stoichiometric amorphous samples obtained at various cooling rates. In the 
inset is shown the relative concentration of Si-Si pairs with respect to Si-N pairs as a function of the cooling rate. 

Fig. 3: (color online) Snapshots of the N(Si 2 N) (A), N(SiN 2 ) (B) and N(N 3 ) (C) units. Circles indicates the 
units containing N-N bonds; N and Si atoms part of these units are bigger and colored in light gray and Orange, 
respectively. 

Fig. 4: Radial Distribution Function (RDF) of the stoichiometric sample obtained from a quenching-from-the-melt 
run at T m = 6000 K and dT/dt = 5 K/ps. In the inset is reported the experimental RDF obtained from neutron 
diffraction, data^ 

Fig. 5: Si™ (top) and N m (bottom) concentration of coordination defect as a function of the cooling rate. 

Fig. 6: Computation vs experimental density at various stoichiometries. Experimental data have been extrapolated 
from data published in Ref. |11) . 

Fig. 7: Total (bottom), Si-Si partial (top left) and N-N partial (top right) pair correlation function of samples at 
various stoichiometries. 

Fig. 8: Concentration of Si™ (top) and N m (bottom) species as a function of the stoichiometry of the sample. 
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TABLE I. List of simulations and corresponding quenching-from-the-melt parameters. Each simulation was followed by a 100 ps 
of thermalization at 300 K and a 10 ps simulation over which the observables have been computed. 



Name 


Melting 
T m [K] tm [ps] 


Cooling 
dT/dt [K/ps] 


5000-100-5 


5000 


100 


5 


5800-100-5 


5800 


100 


5 


(snnn mo i 

UUUU 1UU 1 


6000 


100 


1 


6000-100-5 


6000 


100 


5 


6000-100-10 


6000 


100 


10 


6000-100-50 


6000 


100 


50 


6000-100-100 


6000 


100 


100 


6000-100-500 


6000 


100 


500 


7000-100-5 


7000 


100 


5 



TABLE II. Density used for the NVT simulations as a function of the composition of the sample. 



sample 


X 


density (g/cm 3 ) 


Si64Nl60 


0.400 


2.643 


Sis4Ni4o 


0.600 


2.771 


SiiooNi24 


0.806 


2.890 


Siii2Nn2 


1.000 


2.991 


SiioiNi23 


1.218 


3.091 


Si9eNi2s 


1.333 


3.140 


S193N131 


1.409 


3.169 


Si86Nl38 


1.605 


3.237 


Si89Nl44 


1.800 


3.294 



TABLE III. Concentration and structural data of the most abundant N(Si n N m ) and Si(Si„N m ) units in the melted Si 3 N4 
sample at various T. The structural data reported are the N-N (dNN) and N-Si (dNSi) bond length, and the average bond 
angle (a). These data are computed by time averaging the corresponding instantaneous quantities along the classical and ab 
initio MD simulations. Distances are reported in A and angles in °. 



Classical Simulations 


Ab initio simulations 




5800 K 


6000 K 


7000 K 




5800 K 


6000 K 


7000 K 


N(Sis) & 


ft 


# 


1% 


Si(N 4 ) dfL 


ft 


ft 


ft 


N(Si 2 N) d%° N 


i 






Si(N 3 Si) dfL 

asiSi 
a 








N(Sia) d|k 






w 


Si(N 3 ) dfL 


m 


ft 


ft 


N(SiN) dl° N 








Si(SiN 2 ) dJ° lN 
asiSi 

n 


i 


i 




N(SiN 2 ) dZ° N 

dl^Si 








Si(SiN 4 ) dJ° lN 
asiSi 
a 








N(Si 4 ) d^si 




ft 


ft 


Si(Si 2 N 2 ) dfL 

QSiSi 
a 








N(Si 3 N) 








Si(Si 2 N 3 ) dfL 
a 
















Si(N 2 ) df iN 


ft 




ft. 
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TABLE IV. Concentration and structural data of the most abundant N(Si n N m ) units ( % > 1) in the samples at various 
composition. The structural data reported are the N-N (dNN) and N-Si (dNSi) bond length, and the average bond angle (a). 
These data are computed by time averaging the corresponding instantaneous quantities along the classical and ab initio MD 
simulations. Distances are reported in A and angles in °. 



K - 


= 04 


ic = 


= n r 




X = 


= ns 


X - 


= 1 n 




x = 


= 1 9 


N(Si 3 ) 


A % 

ONSi 

rv 


ft 


N(Si 3 ) 


V 


CO 

ft 


N(Si 3 ) 


A% 


ft 


N(Sia) 


<y, 
j '° 

O-NSi 
IT 




N(Si 3 ) 


^% 
d^si 


ft 


N(bi4j 


A% 

djySi 
a. 


& 




A % 

O-NSi 


w 


JN(bi4j 


A% 




JN(bi4; 


A % 


w 


TYT/O; T\T\ 

IN (012JN J 


A% 






















N(bi 2 N) 


A % 


rfi 








Si 




x = 


= 1 4 




y = 


= 1 fi 


X = 


= 1 8 








N(Si 3 ) 


A % 




N(Sis) 


dl iF i 


ft 


N(Si 3 ) 


A% 


m 


N(Si 2 N) 


^% 

O-NN 
IT 


fa 

119 






N(Si 2 N) 


A % 

Qjvjv 
n 




N(Si 2 N) 


A % 




N(Si 2 N) 


A % 
OiNSi 

O-NN 

rt 




N(Si 3 ) 


A % 

O-NSi 








N(SiN 2 ) 


A % 

% N 




N(SiN 2 ) 


A% 
OiNSi 




N(SiN 2 ) 


A% 
OiNSi 

% N 




N(SiN 2 ) 


A% 
OiNSi 

% N 


i 


















N(N 3 ) 


A% 

O-NN 

rt 




N(N 3 ) 


A% 

rv 


Vis? 










N(Si 3 N) 


A% 
O-NN 

a 















TABLE V. Same data as in Table IV for samples at selected stoichiometries (x = 0.4, Si 3 N4 and x = 1.8) as obtained form 
5 ps long ab initio MD simulations. 



X 


= 0.4 




Si 3 N 4 


x = 


= 1.6 






% 


77.5 




% 


73.5 




% 


38.5 


N(Sis) 


a 


1.82 
117 


N(Sis) 


djvsi 
a 


1.79 
116 


N(Si 2 N) 


d-NSi 

djvjv 

Q 


1.79 
1.42 
115 




% 


22.5 




% 


14.5 




% 


37 


N(Si 4 ) 


diVSi 

a 


1.94 
110 


N(Si 2 N) 


djvsi 
djvjv 
a 


1.80 
1.45 
116 


N(Sis) 


djvsi 

Q 


1.77 
116 










% 


1.5 




% 


12.5 








N(SiN 2 ) 


djvsi 
djviv 
a 


1.83 
1.44 
115 


N(SiN 2 ) 


djvsi 
djvjv 

Q 


1.83 
1.41 
114 














N(N 3 ) 


% 

djvjv 
a 


2 
1.42 
111 
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